Demonstration of one-parameter scaling at the Dirac point in graphene 
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We numerically calculate the conductivity a of an undoped graphene sheet (size L) in the limit 
of vanishingly small lattice constant. We demonstrate one-parameter scaling for random impurity 
scattering and determine the scaling function /3(a) — dlna/dlnL. Contrary to a recent prediction, 
the scaling flow has no fixed point (/? > 0) for conductivities up to and beyond the symplectic metal- 
insulator transition. Instead, the data supports an alternative scaling flow for which the conductivity 
at the Dirac point increases logarithmically with sample size in the absence of intervalley scattering 
— without reaching a scale-invariant limit. 

PACS numbers: 73.20.Fz, 73.20.Jc, 73.23.-b, 73.63.Nm 
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Graphene provides a new regime for two-dimensional 
quantum transport [l|, 0, [H , governed by the absence of 
backscattering of Dirac fermions A counterintuitive 
consequence is that adding disorder to a sheet of undoped 
graphene initially increases its conductivity [fj, Q ■ Inter- 
valley scattering at stronger disorder strengths enables 
backscattering [7] , eventually leading to localization and 
to a vanishing conductivity in the thermodynamic limit 
[1, Q . Intervalley scattering becomes less and less impor- 
tant if the disorder is more and more smooth on the scale 
of the lattice constant a. The fundamental question of 
the new quantum transport regime is how the conductiv- 
ity a scales with increasing system size L if intervalley 
scattering is suppressed. 

In usual disordered electronic systems, the hypothesis 
of one-parameter scaling plays a central role in our con- 
ceptual understanding of the metal-insulator transition 
[ToL [111 ]. According to this hypothesis, the logarithmic 
derivative al In a/d In L = 0(a) is a function only of a it- 
self — irrespective of the sample size or degree of dis- 
order. A positive /3-function means that the system scales 
towards a metal with increasing system size, while a nega- 
tive /3-function means that it scales towards an insulator. 
The metal-insulator transition is at j3 = 0, (3' > 0. In a 
two-dimensional system with symplectic symmetry, such 
as graphene, one would expect a monotonically increas- 
ing /3-function with a metal-insulator transition at (l3j 
as ~ 1.4 (see Fig. [TJ green dashed curve). 

Recent papers have argued that graphene might devi- 
ate in an interesting way from this simple expectation. 
Nomura and MacDonald [l4| have emphasized that the 
very existence of a /3-function in undoped graphene is 
not obvious, in view of the diverging Fermi wave length 
at the Dirac point. Assuming that one-parameter scal- 
ing does hold, Ostrovsky, Gornyi, and Mirlin [l5[ have 
proposed the scaling flow of Fig. Q] (black solid curve). 
Their /3-function implies that a approaches a universal, 
scale invariant value a* in the large- L limit, being the 
hypothetical quantum critical point of a certain field tile- 
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FIG. I: Two scenarios for the scaling of the conductivity a 
with sample size L at the Dirac point in the absence of inter- 
valley scattering. The black solid curve with two fixed points 
is proposed in Ref. []1|, the green dotted curve without a 
fixed point is an alternative scaling supported by the numeri- 
cal data presented in this work. For comparison, we include as 
a red dashed curve the scaling flow in the symplectic symme- 
try class, which has a metal-insulator transition at as ~ 1-4 

m. 



ory. This field theory differs from the symplectic sigma 
model by a topological term [l5l . [HJ . The quantum crit- 
ical point could not be derived from the weak-coupling 
theory of Ref. [l5| , but its existence was rather concluded 
from the analogy to the effect of a topological term in the 
field theory of the quantum Hall effect [ill, S3] • The pre- 
cise value of a* is therefore unknown, but it is well con- 
strained [l3 : From below by the ballistic limit <tq = 1/n 
[HI) [3 H3] and from above by the unstable fixed point 
a s « 1.4. 

In this work we present a numerical test firstly, of the 
existence of one-parameter scaling, and secondly of the 
scaling prediction of Ref. [l5| against an alternative scal- 
ing flow, a positive (3 without a fixed point (green dotted 
curve in Fig. [T|). For such a test it is crucial to avoid 
the finite-a effects of intervalley scattering that might 
drive the system to an insulator before it can reach the 
predicted scale invariant regime. We accomplish this by 
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starting from the Dirac equation, being the a — > limit 
of the tight-binding model on a honeycomb lattice. We 
have developed an efficient transfer operator method to 
solve this equation, which we describe before proceeding 
to the results. 

The single- valley Dirac Hamiltonian reads 

H = vp- a + V{x) + U{x,y). (1) 

The vector of Pauli matrices er acts on the sublattice 
index of the spinor p = —ihd/dr is the momentum 
operator, and v is the velocity of the massless excitations. 
The disorder potential U (r) varies randomly in the strip 
< x < L, < y < W (with zero average, (U) 
I ) : . This disordered strip is connected to highly doped 
ballistic leads, according to the doping profile V(x) = 
for < x < L, V(x) — > — oo for x < and x > L. We 
set the Fermi energy at zero (the Dirac point), so that 
the disordered strip is undoped. The disorder strength is 
quantified by the correlator 

Following Refs. 0, HH, we work with a transfer op- 
erator representation of the Dirac equation H^f = 
at zero energy. We discretize x at the N points 
X\, x^i ■ ■ ■ xn and represent the impurity potential by 
U(r) — J2 n U n (y)6(x — x n ). Upon multiplication by ia x 
the Dirac equation in the interval < x < L takes the 
form 

d 

hv—^ x (y) = [vp y a z - ia x ^2U 7l (y)S(x - x n )]H x (y). 

n 

(3) 

The transfer operator A4, defined by ^ l = -M^q, is 
given by the operator product 

M = Vl,x n ICn'P Xn ^ n _ 1 ■ ■ ■K,2 r P X2:Xl JCi'P Xl fl 1 (4) 
V XjX > = exp[(l/h)(x- x')p y a z ], (5) 
JC n = exp[—(i/hv)U n a x ]. (6) 

The operator V gives the decay of evanescent waves be- 
tween two scattering events, described by the operators 
K n . For later use we note the current conservation rela- 
tion 

M- 1 = a x M^a x . (7) 

We assume periodic boundary conditions in the y- 
direction, so that we can represent the operators in the 
basis 

^ = -=e^\±), q k = W , fc = 0,±l,±2.... (8) 

The spinors |+) = ^ 1,2 (\), R = 2~ 1 / 2 (_ 1 1 ) are eigen- 
vectors of a x . In this basis, (p y )kk' — hqk$kk' is a diago- 
nal operator, while (U n )kk> = W^ 1 J dy U n (y) exp[i(q fc / - 
Qk)y] is nondiagonal. We work with finite-dimensional 



transfer matrices by truncating the transverse momenta 
q k at \k\ = M. 

The transmission and reflection matrices t, r are de- 
termined as in Ref. by matching the amplitudes of 
incoming, reflected, and transmitted modes in the heav- 
ily doped graphene leads to states in the undoped strip 
at x = and x — L. This leads to the set of linear 
equations 

i s kk^4(y) + rkk^kiv)} = My), M 

k 

Y, t ^'rt(y) = *L(y) = M* (y)- (9b) 

Using the current conservation relation ([7]) we can solve 
Eq. ^ for the transmission matrix, 

(l + r) = " Mt (t) = (+\ M *\+)- ( 10 ) 

The transmission matrix determines the conductance 
G = (4e 2 //i) Trttf , and hence the dimcnsionlcss conduc- 
tivity a = (h/Ae 2 )(L/W)G. The average conductivity 
(a) is obtained by sampling some 10 2 — 10 3 realizations 
of the impurity potential. 

Because the transfer matrix V has both exponentially 
small and exponentially large eigenvalues, the matrix 
multiplication ^ is numerically unstable. As in Ref. [22l |. 
we stabilize the product of transfer matrices by trans- 
forming it into a composition of unitary scattering ma- 
trices, involving only eigenvalues of unit absolute value. 

We model the disorder potential U (r) = 5Z n=1 j n S(x— 
x n )5(y — y n ) by a collection of N isolated impurities dis- 
tributed uniformly over a strip 0<x<L,0<y<W. 
(An alternative model of a continuous Gaussian random 
potential is discussed at the end of the paper.) The 
strengths j n of the scatterers are uniform in the inter- 
val [—70,70]- The number TV sets the average separation 
d = (WL/N) 1 / 2 of the scatterers. The cut-off |fc| < M 
imposed on the transverse momenta qk limits the spatial 
resolution £ = W/(2M + 1) of plane waves oc e iqky±q " x at 
the Dirac point. The resulting finite correlation lengths 
of the scattering potential in the x- and y-directions scale 
with £, but they are not determined more precisely. The 
disorder strength ([2]) evaluates to K = ^j 2 (hvd)~ 2 , in- 
dependent of the correlation lengths. We scale towards 
an infinite system by increasing M oc L at fixed disorder 
strength K , scattering range and aspect ratio W/L. 

This completes the description of our numerical 
method. We now turn to the results. In Fig. [2] we first 
show the dependence of the average conductivity on Kq 
for a fixed system size. As in the tight-binding model 
of Ref. @ , disorder increases the conductivity above the 
ballistic value. This impurity assisted tunneling |5j satu- 
rates in an oscillatory fashion for Kq 3> I (unitary limit 
[H Hi). In the tight-binding model [| the initial in- 
crease of a was followed by a rapid decay of the conduc- 
tivity for Kq > 1, presumably due to Anderson localiza- 
tion. The present model avoids localization by eliminat- 
ing intervalley scattering from the outset. 
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FIG. 2: Disorder strength dependence of the average conduc- 
tivity for a fixed system size ( W — 4 L = 40 d) and four values 
of the scattering range. 



The system size dependence of the average conductiv- 
ity is shown in Fig. [3J for various combinations of disorder 
strength and scattering range. We take W/L sufficiently 
large that we have reached an aspect-ratio independent 
scaling flow and L/d large enough that the momentum 
cut-off M > 25. The top panel shows the data sets as 
a function of L/d. The increase of a with L is approx- 
imately logarithmic, (a) — constant + 0.25 In L, much 
slower than the \fL increase obtained in Ref. in the 
absence of mode mixing. 

If one-parameter scaling holds, then it should be pos- 
sible to rescale the length L* = f(Ko,£/d)L such that 
the data sets collapse onto a single smooth curve when 
plotted as a function of L*/d. (The function / = d/l* de- 
termines the effective mean free path I*, so that L* jd = 
L/l* .) The bottom panel in Fig. [^demonstrates that, in- 
deed, this data collapse occurs. The resulting /3-function 
slotted in the inset. Starting from the ballistic limit 
at <7o — l/ 71 "; the /3-function first rises until a « 0.6, 
and then decays to zero without becoming negative. For 
a > as ~ 1-4 the decay oc 1/a is as expected for a 
diffusive system in the symplectic symmetry class. The 
positive /3-function in the interval ((Jo,<js) precludes the 
flow towards a scale-invariant conductivity predicted in 
Ref. 

The model of isolated impurities considered so far is 
used in much of the theoretical literature, whereas exper- 
imentally a continuous random potential is more realistic 
[LH ]. We have therefore also performed numerical simu- 
lations for a random potential landscape with Gaussian 
correlations 12511 . 



{ U{r)U{r')) = K ^%e-^ 'l^ 2 . (11) 
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The discrete points x\, xi . . . xn in the operator product 
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FIG. 3: System size dependence of the average conductivity, 
for W/L = 4 (black and green solid symbols) and W/L — 1.5 
(all other symbols) and various combinations of Kq and £/d. 
The top panel shows the raw data. In the bottom panel the 
data sets have been given a horizontal offset, to demonstrate 
the existence of one-parameter scaling. The inset shows the 
resulting /3-function. 



(|4]) are taken equidistant with spacing 8x = L/N, and 



rx n +5x/2 

U n (y) = / dxU(x,y). 

Jx n -Sx/2 



(12) 



We take M, N, and W/L large enough that the resulting 
conductivity no longer depends on these parameters. We 
then scale towards larger system sizes by increasing L/£ 
and W/l; at fixed K . No saturation of a with increasing 
Kq is observed for the continuous random potential (as 
expected, since the unitary limit is specific for isolated 
scatterers [H, HJ). Fig. [4] shows the size dependence 
of the conductivity — both the raw data as a function 
of L (inset) as well as the rescaled data as a function 
of L* = g(K )L. Single-parameter scaling applies for 
L > 5£, where (a) = constant + 0.32 In L. The prefactor 
of the logarithm is about 25% larger than in the model of 
isolated impurities (Fig. [3]), which is within the numerical 
uncertainty. 
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FIG. 4: System size dependence of the average conductivity 
in the continuous potential model, for several values of Ko- 
The inset shows the raw data, while the data sets in the main 
plot have a horizontal offset to demonstrate one-parameter 
scaling when L > 5£. 



In conclusion, we have demonstrated that the central 
hypothesis of the scaling theory of quantum transport, 



the existence of one-parameter scaling, holds in graphene. 
The scaling flow which we find (green dotted curve in Fig. 
[IJ is qualitatively different both from what would be ex- 
pected for conventional electronic systems (red dashed 
curve) and also from what has been predicted [24[ for 
graphene (black solid curve). Our scaling flow has no 
fixed point, meaning that the conductivity of undoped 
graphene keeps increasing with increasing disorder in the 
absence of intervalley scattering. The fundamental ques- 
tion "what is the limiting conductivity Coo of an infinitely 
large undoped carbon monolayer" has therefore three dif- 
ferent answers: doo = 1/-7T in the absence of any disorder 
.18, 19]. doc = oo with disorder that does not mix the val- 
leys (this work), and (Too = with intervalley scattering 

||. 
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Note added: Since submission of this manuscript, simi- 
lar conclusions have been reported by Nomura, Koshino, 
and Ryu 
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